function [gbest,gbestval,fitcount] = SGSA_func(fhd,Dimension,Particle_Number,Max_Gen,lowerbound,upperbound,pos,bool_Plot,varargin)

G0 = 2.8;
epsilon = 1e-3;
rand('state',sum(100*clock));
me = Max_Gen;
ps = Particle_Number;
D = Dimension;

[pub, plb, vub, vlb] = SetBound(D,ps,lowerbound,upperbound);
% pos = VRmin+(VRmax-VRmin).*rand(ps,D);

fitcount = ps;
vel = zeros(ps,D); % initialize the velocity of the particles
acc = zeros(ps,D); % initialize the acceleration of the particles

fitness = feval(fhd,pos',varargin{:});
sort_fitness = [fitness' (1:ps)'];
sort_fitness = sortrows(sort_fitness, 1);
fitbest = sort_fitness(1,1);
fitworst = sort_fitness(ps,1);
gbestindex = sort_fitness(1,2);
gbest = pos(gbestindex,:);
gbestval = fitbest;
if(bool_Plot)
    figure(2)
    haxes = plot( 0 , 0 );
    XArray = [1];
    YArray = [gbestval];
    title('SGSA');
end

for iter = 2:me
    if(fitbest == fitworst)
        mass = repmat(1/ps,ps,1);
    else
        mass = (fitness'-fitworst)/(fitbest-fitworst);
    end
    mass = mass / sum(mass);

    G = G0 * exp(1 - (iter / me));
    numKbest = ceil(ps * (1 - (iter /me)));
    numKbest = (numKbest>1)*numKbest + (numKbest<=1)*1;
    for n = 1 : ps
        acc(n,:) = zeros(1,D);
        for k = 1 : numKbest
            force = zeros(1,D);
            indexKbest = sort_fitness(k, 2);
            force(1,:) = pos(indexKbest,:) - pos(n,:);
            R = norm(force(1,:));
            force(1,:) = force(1,:) * rand() * mass(indexKbest) / (R + epsilon);
            acc(n,:) = acc(n,:) + force(1,:);
        end
    end
    acc = G * acc;
    vel = vel.*rand(ps,D) + acc;
    vel = ((vel>=vlb)&(vel<=vub)).*vel+(vel<vlb).*vlb+(vel>vub).*vub;
    pos = pos + vel;
%     pos = ((pos>=plb)&(pos<=pub)).*pos+(pos<plb).*plb+(pos>pub).*pub;
    pos=((pos>=plb)&(pos<=pub)).*pos+(pos<plb).*(plb+0.25.*(pub-plb).*rand(ps,D))+...
        (pos>pub).*(pub-0.25.*(pub-plb).*rand(ps,D));
    
    fitness = feval(fhd,pos',varargin{:});
    sort_fitness = [fitness' (1:ps)'];
    sort_fitness = sortrows(sort_fitness, 1);
    fitbest = sort_fitness(1,1);
    fitworst = sort_fitness(ps,1);
    gbestindex = sort_fitness(1,2);
    gbest = pos(gbestindex,:);
    gbestval = fitbest;
    
    if(bool_Plot & mod(iter,10)== 0)
        XArray = [ XArray iter]; 
        YArray = [ YArray gbestval];
        set( haxes , 'XData' , XArray , 'YData' , YArray );
        drawnow
    end
end

end
